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Abstract.  The  convergence  behavior  of  Bird’s  Direct  Simulation  Monte  Carlo  (DSMC)  method  is  explored  for  near¬ 
continuum,  one-dimensional  Fourier  flow.  An  argon-like,  hard-sphere  gas  is  confined  between  two  parallel,  fully 
accommodating,  motionless  walls  of  unequal  temperature.  The  simulations  are  performed  using  a  one-dimensional  code 
based  on  Bird’s  DSMC  prescription.  The  convergence  metric  is  taken  as  the  ratio  of  the  DSMC-calculated  bulk  thermal 
conductivity  to  the  infinite-approximation,  Chapman-Enskog  (CE)  theoretical  value.  The  convergence  rate  is  monitored 
with  respect  to  three  parameters:  time  step,  cell  size,  and  number  of  computational  molecules  per  cell.  For  a  sufficiently 
large  number  of  computational  molecules,  the  DSMC  discretization  error  is  observed  to  decrease  quadratically  with  time 
step  and  cell  size,  in  good  agreement  with  previous  predictions  based  on  Green-Kubo  theory.  When  all  three  parameters  are 
finite,  the  observed  convergence  behavior  is  complex.  As  numerical  errors  are  systematically  reduced,  the  DSMC- 
calculated  conductivity  is  shown  to  approach  the  theoretical  CE  value  to  within  0.02%,  showing  that  the  present 
implementation  of  Bird’s  DSMC  algorithm  provides  an  excellent  representation  of  continuum,  CE  heat  conduction. 


INTRODUCTION 

The  Direct  Simulation  Monte  Carlo  (DSMC)  method  employs  a  particle -based,  stochastic  algorithm  to  solve  the 
Boltzmann  equation  by  approximating  the  continuous  molecular  velocity  distribution  function  with  a  discrete  set  of 
computational  molecules  [1].  Historically,  DSMC  has  been  very  successful  in  the  study  of  rarefied,  high-speed  flows 
typical  of  aerospace  engineering  [1,2].  The  validity  of  the  method  for  this  class  of  problems  has  been  firmly  established 
by  comparisons  to  experimental  data  [3]  and  molecular-dynamics  calculations  [4].  More  recently,  the  DSMC  method 
has  become  increasingly  popular  as  a  means  for  exploring  ambient-pressure  flows  in  nano-  and  micro-scale 
geometries,  such  as  are  found  in  thin-film  bearings  [5],  oscillating  microbeams  [6],  and  microchannels  [7,8].  Often, 
these  small-scale  systems  involve  near-continuum  and/or  low-speed  flows  that  are  challenging  to  simulate  with 
stochastic  algorithms.  For  example,  one  DSMC  study  of  ambient-pressure,  low-speed  (~1  m/s)  micro-Couette  flow 
suggests  that  statistical  errors  become  large  enough  to  preclude  accurate  simulations  [8].  Near-continuum  DSMC 
studies  of  the  motionless-gas  Fourier  problem  have  been  more  successful.  For  example,  several  studies  report  DSMC 
simulations  that  give  bulk  thermal  conductivities  that  agree  with  Chapman-Enskog  (CE)  theoretical  predictions  to 
within  -10%  [9,10].  These  authors  note  that  some  of  the  observed  deviations  may  arise  from  the  finite  nature  of  the 
computational  domain.  Recently,  Gallis  et  al.  [11]  performed  very  careful  DSMC  simulations  of  the  Fourier  problem 
and  found  excellent  agreement  (-0.2%)  between  their  calculated  bulk  thermal  conductivity  and  the  CE  infinite- 
approximation  result.  They  also  showed  that  the  DSMC-calculated  velocity  distribution  function,  as  measured  by  the 
Sonine-polynomial  coefficients  Ujt,  was  in  excellent  agreement  with  CE  theoretical  predictions  [12],  Thus,  at  least  for 
the  Fourier  problem,  it  appears  that  careful  applications  of  the  DSMC  method  can  provide  solutions  that  are  in  good 
agreement  with  the  CE  approximation  to  the  Boltzmann  equation. 

One  drawback  to  achieving  high-accuracy  solutions  with  the  DSMC  method  is  the  high  computational  burden. 
Consequently,  a  clear  understanding  of  how  to  achieve  a  specified  level  of  numerical  accuracy  with  a  minimum  of 
computational  effort  is  needed.  In  response,  this  paper  provides  a  convergence  study  of  the  DSMC  method  for  the  one¬ 
dimensional  Fourier  heat-flow  problem.  The  convergence  metric  is  taken  as  the  ratio  of  the  DSMC-calculated  bulk 
thermal  conductivity,  KDSMC  to  the  infinite-approximation,  CE  theoretical  value,  K  [11,12].  To  ensure  that  the  CE 
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limit  is  achieved,  DSMC  simulations  are  performed  for  a  small  system  Knudsen  number  (-0.02)  and  with  a  modest 
temperature  gradient  (100  K/mm).  These  conditions  lead  to  solutions  for  which  the  normal  solution  in  the  central 
region  of  the  domain  can  be  clearly  differentiated  from  the  Knudsen  layers  near  each  wall.  Four  parameters  are  known 
to  limit  the  numerical  accuracy  of  the  DSMC  method:  the  number  of  independent  samples  per  cell,  Mc\  the  number  of 
simulators  (particles)  per  cell,  Nc\  the  time  step.  At;  and  the  cell  size.  Ax.  The  present  calculations  employ  large  enough 
samples  such  that  statistical  errors  are  reduced  to  levels  that  are  negligible  compared  to  the  contributions  of  the  other 
three  parameters.  The  remaining,  non-statistical  error  (herein  referred  to  as  the  discretization  error)  is  systematically 
explored  as  a  function  of  Nc,  At,  and  Ax. 


FOURIER  PROBLEM 


The  Fourier  problem  consists  of  one-dimensional  heat  flow  through  a  hard-sphere  gas  confined  between  two 
parallel,  fully  accommodating,  motionless  walls  of  unequal  temperature.  The  wall  at  x  =  0  is  held  at  a  lower 
temperature,  T ,  than  the  temperature,  T H  ,  of  the  wall  at  x  =  L  (Figure  1).  Because  of  its  geometric  simplicity, 
the  Fourier  problem  has  been  the  subject  of  many  studies.  One  common  observation  has  been  that  the  linear 
relationship  between  the  heat  flux,  q ,  and  the  local  temperature  gradient  (Fourier’s  law),  given  in  one  dimension  by 


dT 

q  =  -UT)-. 


(1) 


remains  valid  even  under  strongly  nonequilibrium  (high  heat  flux)  conditions.  The  validity  of  Equation  1  for  high  heat 
fluxes  has  been  predicted  theoretically  for  gases  undergoing  Maxwell  [13]  and  BGK  [14]  interactions,  and  has  also 
been  demonstrated  for  hard-sphere  gases  using  both  molecular  dynamics  [15,16]  and  DSMC  [9,10,11]  numerical 
simulations.  In  particular,  these  studies  show  that  the  thermal  conductivity  observed  at  high  heat  flux  equals  that 
obtained  from  CE  theory  for  a  nonequilibrium  gas  in  the  continuum  limit  of  small  heat  flux  and  Knudsen  number.  CE 
theory  for  a  hard-sphere  gas  relates  the  thermal  conductivity  to  a  reference  viscosity,  | xre^ ,  as  [12,1 1]: 
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where  Krej  and  (i rej  are  evaluated  at  the  reference  temperature  T  rej ,  kB  =  1.380658x10  -  J/K  is  the  Boltzmann 
constant,  K^/K^  and  are  known  CE  values,  and  m  is  the  molecular  mass.  For  a  hard-sphere  gas,  Gallis  et 

al.  [11]  used  CE  theory  [12]  to  calculate  the  infinite-order  values  /K  j  =  1.025218  and  p^/pi  =  1.016034. 

Gallis  et  al.  [11]  performed  DSMC  simulations  that  showed  that  Fourier’s  law  remains  valid  at  high  heat  fluxes 
despite  the  fact  that  the  underlying  molecular  velocity  distribution  deviates  from  CE  theory,  as  measured  by  the 
departure  of  the  DSMC-calculated  Sonine -polynomial  coefficients  from  the  expected  CE  values.  Their  work  showed 
that  the  DSMC  method  reproduces  the  CE  velocity  distribution  only  in  the  limit  of  small  local  Knudsen  numbers 
(<  0.01);  at  higher  values,  the  Sonine -polynomial  coefficients  rapidly  diverge  from  the  CE  prediction. 


DISCRETIZATION  ERRORS  IN  THE  DSMC  METHOD 

Wagner  [17]  presented  a  rigorous  mathematical  proof  that  numerically  resolved  (i.e.,  in  the  limit  of  infinite  Nc) 
DSMC  simulations  provide  solutions  of  the  Boltzmann  equation.  This  result  suggests  that  calculations  performed  with 
a  careful  implementation  of  the  DSMC  method  should  agree  with  theoretical  solutions  of  the  Boltzmann  equation. 
Four  parameters  control  numerical  error  in  the  DSMC  method  [1]:  the  sample  size,  Mc,  the  number  of  simulators  per 
cell,  Nc ,  the  time  step.  At,  and  the  cell  size.  Ax.  Bird  [1]  observed  that  statistical  fluctuations  decline  with  the  square 
root  of  the  sample  size  and  could  be  reduced  (in  principle)  to  any  desired  level  by  continuing  or  repeating  the 
simulation.  Fladjiconstantinou  et  al.  [18]  recently  presented  closed-form  expressions  that  relate  the  statistical  error  to 
the  sample  size  in  molecular  simulation  algorithms.  In  the  present  study,  sufficiently  large  samples  are  taken  so  that 
statistical  errors  are  reduced  to  negligible  levels. 
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FIGURE  1.  Temperature  profile  for  the  Fourier  geometry  (inset).  FIGURE  2.  Conductivity-ratio  profile. 


Bird  distinguished  statistical  error  from  the  discretization  errors  introduced  by  the  three  remaining  parameters  (Nc, 
At,  and  Ax).  Early  practitioners  of  the  DSMC  method  relied  on  rule-of-thumb  guidelines  to  limit  discretization  errors. 
For  example.  Bird  [1]  recommended  that  ~20  or  more  simulators  per  cell  be  used,  that  cell  dimensions  should  be  kept 
below  —1/3  of  the  local  mean  free  path,  and  that  the  time  step  should  be  kept  below  -1/4  of  the  local  mean  collision 
time.  Later  studies  have  sought  to  quantify  the  discretization  errors  associated  with  the  various  numerical 
approximations.  Chen  and  Boyd  [19]  analyzed  the  numerical  error  inherent  in  the  DSMC  method  in  terms  of  the 
sample  size  and  the  number  of  simulators  per  cell.  Their  calculations  show  that  the  error  contains  two  contributions,  a 
statistical  error  proportional  to  I  /  JM  r  and  a  discretization  error  proportional  to  1  /Nc.  Thus,  in  the  limit  of  vanishing 
time  step  and  cell  size  (and  infinite  sample  size),  the  discretization  error  in  the  ratio  of  the  DSMC-calculated  thermal 
conductivity  to  the  CE  theoretical  value  is  expected  to  be  of  the  form: 
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where  A  is  a  constant  that  has  not  been  previously  reported  in  the  literature. 

Alexander  et  al.  [20]  applied  Green-Kubo  theory  to  derive  explicit  expressions  for  the  dependence  of  hard-sphere 
gas  transport  coefficients  on  cell  size.  Their  expression  for  the  thermal-conductivity  ratio,  in  the  limit  of  vanishing  time 
step  and  infinite  simulators  per  cell,  is: 

,.  kdsmc  ,  32 
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where  Ax  is  width  of  the  cell  in  the  x-direction,  X  -  ( J2nd^ej-n)~l  is  the  hard-sphere  gas  mean  free  path,  d rgj-  is  the 
molecular  diameter,  and  n  is  the  number  density.  DSMC  calculations  [20]  of  the  momentum  flux  for  the  Couette-flow 
problem  confirmed  the  predicted  0(  Ax2 )  truncation  error  and  agreed  within  numerical  uncertainty  with  the  Green- 
Kubo  predictions  although  only  two  data  points  were  reported  in  the  critical  range  of  interest  Ax/X  <  1  . 

Hadjiconstantinou  [21]  applied  Green-Kubo  theory  to  derive  explicit  expressions  for  the  dependence  of  hard- 
sphere  gas  transport  coefficients  on  time  step.  His  expression  for  the  thermal-conductivity  ratio,  in  the  limit  of 
vanishing  cell  size  and  infinite  simulators  per  cell,  is: 
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where  At  is  the  time  step,  t()  =  'k/cQ  is  the  hard-sphere  mean  collision  time  and  ca  =  J2k ^T/m  is  the  most 
probable  molecular  speed.  DSMC  simulations  for  a  variety  of  one-dimensional  flows  [22]  confirmed  the  predicted 
0(At2 )  truncation  error  although  only  a  few  data  points  were  reported  in  the  critical  range  of  interest  A t/t  <  1  . 


DSMC  SIMULATIONS 


The  convergence  behavior  of  the  DSMC  method  is  determined  for  one-dimensional  Fourier  heat  flow  of  a  hard- 
sphere  gas  between  parallel  walls  of  unequal  temperature  (see  Figure  1  inset).  The  cold  and  hot  walls  are  held  at  223.15 
and  323. 15  K,  respectively,  and  are  separated  by  a  distance  of  L  =  1  mm  .  All  DSMC  results  below  are  acquired  using 
a  modified  version  of  Bird’s  one-dimensional  code  DSMC1  [1,11].  The  main  modification  is  the  recasting  of  the  code 
into  double  precision.  In  keeping  with  the  spirit  of  high  accuracy,  round-off  tolerances  in  several  routines  were 
eliminated.  The  strategy  of  sampling  molecules  both  before  and  after  the  collision  stage  (move-sample-collide-sample, 
MSCS)  was  used  [11],  A  hard-sphere  gas  is  modeled  using  argon-like  properties  (m  =  6.63x10"  kg, 
\Xref  =  2.117x10  5  Pa  •  s  at  T rej  =  273.15  K)  with  the  VSS  collision  model  (to  =  0.5  and  a  =  1 ).  The  molecular 
diameter  is  calculated  using  the  Bird  strategy  [1]  as  modified  by  Gallis  et  al.  [11]  to  account  for  higher-order  CE 
corrections;  the  present  value  for  argon-like  hard  spheres  is  dref  =  3.658x10  m  .  In  all  cases,  the  gas  is  originally 
motionless  at  temperature  T  rep  and  pressure  p rej  =  266.44  Pa ,  and  flow  transients  are  allowed  to  decay  for  7  ms 
before  the  steady  state  phase  is  initiated.  The  final  DSMC-calculated  pressure,  P dsmc  ~  266.9  Pa ,  differs  slightly 
from  the  initial  pressure.  The  final  pressure  and  the  reference  temperature  are  used  to  specify  the  mean  free  path  and 
collision  time:  A.(7’rgy,  =  2.40x10  m  and  tg{T  rep  P  dsmA  =  7 1  ns  .  The  domain  is  ~42  mean  free  paths 

across,  indicating  that  the  Knudsen  layers  at  the  walls  occupy  only  a  small  fraction  of  the  domain. 

Time  step,  cell  size,  and  the  number  of  simulators  per  cell  are  systematically  varied.  The  dimensionless  time  step 
was  selected  from  the  range  0.05  <  A  t/t  <  1 ,  corresponding  to  physical  time  steps  between  3.5  and  70  ns.  The 
dimensionless  cell  size  was  selected  from  the  range  0.05  <  Ax/'k.  <  0.85  ,  corresponding  to  800  to  50  cells, 
respectively,  across  the  domain.  The  number  of  simulators  used  in  each  cell  was  selected  from  the  range  7  <  N  <  240  . 
Each  simulation  is  run  for  ~48  hours  on  16  nodes  of  an  IBM  Linux  cluster,  with  dual  1.2-GHz  P3  processors  per  node. 
The  parallel  implementation  of  the  code  speeds  the  calculations  and  uses  ensemble  averaging  to  provide  95%- 
confidence  intervals  for  all  calculated  values.  Runs  are  continued  until  this  uncertainty  becomes  small  compared  to  the 
discretization  errors,  typically  requiring  O(1010)  samples  per  cell.  KDSMC/K  is  calculated  for  each  cell  from 


kdsmc 

K 


' Tref  \  \g\  VdTY1 
KrefjyTV2Adx) 


(6) 


where  q,  T,  and  dT/dx  are  determined  by  DSMC.  The  DSMC-calculated  wall  heat  flux  is  used  for  q  in  Equation  6. 
Additional  numerical  error  is  introduced  if  the  local,  cell-based,  volume-averaged  heat-flux  moment  is  used  instead  of 
the  wall  value.  Note  that  the  present  MSCS  sampling  strategy  improves  the  accuracy  of  volume-averaged  flux 
moments  compared  to  the  standard  strategy  (move-collide-sample,  MCS)  but  has  no  effect  on  wall-based  values.  Thus, 
the  present  convergence  study  applies  equally  well  to  either  the  MCS  or  MSCS  sampling  strategies. 

The  local  cell-based  value  for  temperature  is  used.  A  profile  of  a  typical  DSMC-calculated  temperature  profile  (200 
cells,  7-ns  time  step,  30  simulators  per  cell  )  is  shown  in  Figure  1,  demonstrating  the  expected  temperature  jumps  near 
the  walls  and  a  nearly  linear  profile  in  the  central  region  of  the  domain.  The  temperature  gradient  is  determined  from 
the  temperature  profile  with  a  central-difference  scheme.  A  conductivity-ratio  profile  is  shown  in  Figure  2  which  has 
been  calculated  from  Equations  2  and  6  using  DSMC  simulation  results  for  the  same  parameter  set  as  in  Figure  1 .  The 
ratio  is  observed  to  lie  very  close  to  unity  throughout  the  interior,  with  Knudsen  layers  evident  near  the  walls.  Some 
statistical  scatter  in  the  points  is  observed,  which  is  amplified  by  the  need  to  differentiate  the  temperature  profile.  A 
single  measure  of  the  conductivity  ratio  is  found  for  each  calculation  by  averaging  the  cell-based  conductivity  ratios 
over  the  central  40%  of  the  domain  (well  outside  the  Knudsen  layers);  this  averaging  procedure  greatly  reduces  the 
statistical  scatter.  The  error  bars  for  each  data  point  plotted  in  Figures  3  and  4  are  the  95%  confidence  interval  for  each 
center-averaged  value  of  the  conductivity  ratio. 
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FIGURE  3.  Conductivity-ratio  dependence  on  dimensionless 
time  step  for  various  Nc  (200  cells). 


FIGURE  4.  Conductivity-ratio  dependence  on 
dimensionless  cell  width  for  various  Nc  {At  =  7  ns). 


The  convergence  behavior  of  the  averaged  conductivity  ratios  is  explored  systematically  with  over  200  simulations 
corresponding  to  the  ranges  given  above  for  cell  size,  time  step,  and  number  of  simulators.  One  set  of  calculations  is 
shown  in  Figure  3,  where  the  conductivity  ratio  is  plotted  against  the  square  of  the  dimensionless  time  step  for 
15  <  N  <  240  and  for  200  cells  ( Ax/A,  =  0.209  ).  The  theoretical  prediction  of  Equation  5  is  also  plotted  and  appears 
as  a  straight  line.  For  the  largest  number  of  simulators  considered,  N  =  240 ,  the  DSMC  calculations  lie  parallel  to, 
but  displaced  slightly  upward  from,  the  theoretical  line.  This  displacement  results  from  cell-size  and  finite-simulator 
discretization  errors  which  have  not  been  eliminated.  Although  quadratic  convergence  in  time  step  is  approached  with 
increasing  simulators,  the  convergence  behavior  for  fewer  than  -60  simulators  per  cell  is  clearly  not  quadratic. 
Extrapolation  of  the  entire  convergence  data  set  to  vanishing  cell  size  and  infinite  Nc  gives  an  estimate  for  the  time- 
step  coefficient,  0.0288,  that  is  in  good  agreement  with  the  theoretical  value,  0.03018. 

The  convergence  behavior  of  the  averaged  conductivity  ratio  with  cell  size  is  demonstrated  in  Figure  4;  these 
calculations  were  all  performed  with  a  fixed  time  step  A t/t  =0.1  and  for  15  <  Nc  <  240  .  The  theoretical  prediction 
of  Equation  4  is  also  plotted  and  appears  as  a  straight  line.  Quadratic  convergence  in  cell  size  is  approached  for  the 
largest  Nc  values  tested,  although  the  DSMC  calculations  are  not  as  parallel  to  the  theoretical  result  as  for  the  time- 
step  results.  Extrapolation  of  the  entire  convergence  data  set  to  vanishing  time  step  and  infinite  Nc  gives  an  estimate 
for  the  quadratic  time-step  coefficient,  0.0401,  that  is  in  fair  agreement  with  the  theoretical  value,  0.04527.  Again,  the 
convergence  behavior  with  cell  size  for  fewer  than  -60  simulators  per  cell  is  clearly  not  quadratic. 

A  Taylor  series  expansion  in  terms  of  the  three  parameters  was  used  to  perform  a  least-squares  fit  of  the  entire 
convergence  data  set.  The  result  of  the  fitting,  shown  in  both  Figures  3  and  4  as  dashed  lines,  agrees  extremely  well 
with  the  DSMC  calculations.  The  statistically-significant,  leading  terms  of  the  fitting  are  found  to  be: 
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The  leading  term  differs  from  unity  by  less  than  0.02%,  showing  that  the  present  DSMC  implementation  is  in 
remarkable  agreement  with  CE  theory.  The  coefficients  of  the  l/Nc  and  At  (Ax)2  terms  have  not  been  previously 
reported;  interestingly,  both  reduce  the  conductivity  ratio.  Higher-order  cross  terms  are  also  evident. 


CONCLUSIONS 

A  comprehensive  study  of  the  numerical  convergence  behavior  of  the  DSMC  method  has  been  completed  for  the 
Fourier  problem.  The  study  shows  that  the  present  implementation  of  the  DSMC  method  converges  to  within  0.02% 
of  the  CE  thermal  conductivity  in  the  limit  of  vanishing  cell  size  and  time  step  with  infinite  simulators.  Quadratic 
convergence  in  time  step  and  cell  width  is  observed  for  N  — »  °° ,  in  good  agreement  with  Green-Kubo  theory. 
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